%%TORQSTOCK_BUILD.m
% This program creates TORQSTOCK, which will be used in Stata

% This program is part of the data and code used in:
% Kenneth R. Ahern, "Do Common Stocks Have Perfect Substitutes? 
% Product Market Competition and the Elasticity of Demand for Stocks," 
% The Review of Economics and Statistics, October 2014, 96(4):756-766.

clear all
global CRM CRD CSD

load CQ % original TORQ data
%Get list of symbols from CQ
S1=sortrows(CQ(find(CQ(:,2)==901101),1));
S2=[1;S1(1:end-1)];
SYMBOLS=S1(S1~=S2);
clear S1 S2 CQ
TORQSTOCK = NaN*ones(144,23); % Make placeholder
TORQSTOCK(:,1)=SYMBOLS; %SYMBOLS are the ticker symbols converted to ASCII
clear SYMBOLS
IDS=xlsread('TORQPERMNOGVKEY.xls'); % Data of PERMNOs and GVKEYs
TORQSTOCK(:,2:3)=IDS; % Put PERMNO and GVKEY in TORQSTOCK data
clear IDS

%% Find CRSP variables and add to TORQSTOCK
load FACTORS % Fama-French Daily factors (used to get trading dates)
CRD=csvread(['TORQCRSPDAILYCSV1.csv']); % Stock market data on TORQ firms
clear SDDR1Y DR1 DM1 SHROUT TURNOVER MNDV1Y
%%Get daily dates
LDIND=max(find(FACTORS(:,1)<19901101));
FDIND1=LDIND-250+1;
LD=FACTORS(LDIND,1); % Last date index
FD1=FACTORS(FDIND1,1); % First date index

%%firm loop starts
for i=1:144 %144 firms in TORQ data
    clear SDDR1Y DR1 DM1 SHROUT TURNOVER MNDV1Y
    DR1=CRD(find((CRD(:,1)>=FD1)&(CRD(:,1)<=LD)&(CRD(:,2)==TORQSTOCK(i,2))),[10 9]); %Daily 1 yr Returns/Vol for firm i
    LDR1=length(DR1(:,1));
    if LDR1<250; %Fill in missing dates with NaN
        DM1=NaN*ones(250-LDR1,2);
        DR1=[DM1;DR1];
    end
    if isempty(find(DR1(:,1)<-1))==0;
        DR1(find(DR1(:,1)<-1),:)=NaN*ones(1,2);
    end
    DNM1=sum(isnan(DR1(:,1))==0);  %Number of nonmissing observations
    %%%1 year returns daily
    DRNM1=DR1(find(isnan(DR1(:,1))==0),1);
    SDDR1Y=std(DRNM1);                    %Standard deviation of returns
    %1 year volume daily
    DVNM1=DR1(find(isnan(DR1(:,2))==0),2);
    MNDV1Y=mean(DVNM1);                   % Trading volume
    SHROUT=abs(CRD(max(find((CRD(:,1)>=19891101)&(CRD(:,1)<19901101)&(CRD(:,2)==TORQSTOCK(i,2)))),5)); % Shares outstanding
    TURNOVER=MNDV1Y/SHROUT; %Turnover
    TORQSTOCK(i,4:5)=[SDDR1Y TURNOVER]; % add to TORQSTOCK data
end % end firm loop

%% Find Compustat variables and add to TORQSTOCK
CSD=xlsread('TORQCompustatXLS2.xls'); % Compustat data for TORQ stocks
%%Match on GVKey and for correct date
SIC=NaN(144,1);
for i=1:144 % firm loop
    clear IND1
    IND1=max(find((CSD(:,1)==TORQSTOCK(i,3))&(CSD(:,2)<19901101)));
    if isempty(IND1)==0;
        TORQSTOCK(i,6)=log(CSD(IND1,6));
        SIC(i,1)=CSD(IND1,31); % 4-DIGIT SIC CODE
    end
    clear IND1
end % end firm loop


%*************************************************************
%%ELAS2TORQSTOCK
%% This program is called inside DATA2TORQSTOCK
%% Find stock market elasticity and add to TORQSTOCK
load('TORQELASTICITYEST.mat'); % load elasticity estimates
EDA=NaN*zeros(144,2,2);
for i=1:144
    for b=1:2
        %%Elasticity of demand
        ED=reshape(squeeze(E(i,:,:,b)),63*13,1);  %(Stock:Day:Time) for Buyside
        ND=reshape(squeeze(N(i,:,:,b)),63*13,1);
        COND=reshape(squeeze(CON(i,:,:,b)),63*13,1);
        %%Get daily elasticity by stock (require at least 5 obs [16th percentile])
        for j=1:2
            clear CC ER Q T L l w_l t T2 SE MCC TCC MNOBS

            if j==1;
                CC=ED; %Coefficient estimates over time
            else
                CC=COND;
            end

            CC(ND<5)=NaN; % Set coeffiecients as NaN if less than 5 obs
            ER=CC-nanmean(CC);  % e resid from mean of coefficients
            Q=0;
            T=length(CC);  %total time series length
            L=floor(4*(T/100)^(2/9));  %Lag length according to a Bartlett Kernel and from Newey West 1994
            for l=0:L  %Newey West Autocorrelation Consistent Estimate of Variance
                w_l=1-l/(L+1);
                for t=l+1:T;
                    if (isnan(ER(t))==0)&(isnan(ER(t-l))==0);
                        if l==0 
                            Q=Q+ER(t)^2;
                        else
                            Q=Q+2*w_l*ER(t)*ER(t-l);
                        end
                    end
                end
            end
            if Q>0;
                T2=sum(isnan(CC)==0);
                Q=1/(T2-1)*Q;
                SE=sqrt(Q/T2); %standard error of coefficient (Newey West)
                MCC=nanmean(CC);%average coefficient over time series
                TCC=MCC/SE;
                MNOBS=nanmean(ND(ND>=5));
                EDA(i,j,b)=MCC;
            end
            clear CC ER Q T L l w_l t T2 SE MCC TCC MNOBS
        end  %j elast = 1, cons = 2
    end  %b buy =1, sell = 2
end  %i firms 1...144
ABSEDA=abs(EDA(:,1,1));
TORQSTOCK(:,7)=ABSEDA;

%% Find PIN and add to TORQSTOCK
load PIN8301.txt
PINMR=NaN*ones(144,1);
for i=1:144;
    if isempty(find((PIN8301(:,1)<=1990)&(PIN8301(:,2)==TORQSTOCK(i,2)), 1, 'last' ))==0;
        PINMR(i,1)=PIN8301(find((PIN8301(:,1)<=1990)&(PIN8301(:,2)==TORQSTOCK(i,2)), 1, 'last' ),3); %Most recent PIN thru 1990
    end
end
TORQSTOCK(:,8)=PINMR;

%*************************************************************
%% Find Fama-French Industry 5s and add to data
load SICCODES5MATLAB.txt
FFIND5=zeros(144,1);
for i=1:length(SICCODES5MATLAB);
    if isempty(find((SIC>=SICCODES5MATLAB(i,2))&(SIC<=SICCODES5MATLAB(i,3))))==0;
        FFIND5(find((SIC>=SICCODES5MATLAB(i,2))&(SIC<=SICCODES5MATLAB(i,3))))=SICCODES5MATLAB(i,1);
    end
end
FFIND5(find(FFIND5(:)==0))=5;
FFIND5(find(SIC==1))=NaN;
TORQSTOCK(:,9)=FFIND5;

%% MRR and Lambda
% Unfortunately, the gmm optimization code does not always converge for every observation depending on a random seed.
% Since, different observations don't converge from when I originally calculated these measures, I present the measures from the paper here.
% The matlab .m file called MRRLAMBDA.m presents the code that generated
% the data.  

% MRRRGMM
TORQSTOCK(:,10)=[0.43321
0.77217
0.53638
NaN
0.1423
0.20991
1.0176
0.053796
0.92058
0.076888
0.39822
0.25954
0.22253
0.40163
0.15646
0.62622
0.36427
0.49197
0.32092
0.1419
0.70818
0.55918
0.77706
0.39809
0.54307
0.2191
0.69708
0.19203
0.50712
0.47321
0.73005
-0.22144
0.22656
0.19085
0.5488
0.63718
0.4861
0.13225
NaN
0.63318
0.16646
NaN
NaN
0.27414
NaN
0.40417
0.51465
0.40074
NaN
NaN
0.17206
0.25455
0.14406
0.11301
0.26515
0.14801
NaN
0.034327
0.7609
0.5026
0.26237
0.049775
0.6718
0.41798
0.82718
0.41546
0.06539
0.48203
0.39799
0.55626
0.42105
0.61796
NaN
0.26569
0.18073
0.41324
0.58644
0.39527
1.2286
NaN
NaN
1.0152
0.38317
0.57123
0.57952
0.034968
0.45473
0.5528
1.1543
0.53029
0.25402
0.18393
0.35535
NaN
NaN
0.15087
0.30716
0.37871
0.31743
NaN
0.31126
0.33576
0.079199
0.37407
0.4488
0.46658
0.11782
0.3115
0.44222
0.25971
0.8923
0.43953
0.15217
0.26323
0.46278
NaN
0.42365
NaN
0.23564
0.016711
1.1197
0.29114
0.37485
0.69932
0.55153
0.15869
0.25527
0.46648
0.63467
0.64628
0.39403
0.99673
0.34351
0.21022
NaN
0.27936
0.33235
NaN
0.91413
0.54625
NaN
1.3205
0.47782
0.46971
];

% BSGHLAMBDA
TORQSTOCK(:,11)=[

0.0000121000
0.0000021800
0.0000428000
NaN
0.0000001550
0.0000026500
0.0000069900
0.0000001860
0.0000502000
0.0000005510
0.0000014200
0.0000021000
0.0000003390
0.0000019400
0.0000007300
0.0000037400
0.0000023000
0.0000017800
0.0000010700
0.0000003980
0.0000030400
0.0000014200
0.0000004280
0.0000001930
0.0000102000
0.0000013900
0.0000037200
0.0000004360
0.0000026000
0.0000037300
0.0000021100
0.0000131000
0.0000012300
0.0000013900
0.0000017300
0.0000007130
0.0000023500
0.0000004140
NaN
0.0000071200
0.0000028200
NaN
NaN
-0.0000001200
0.0000004930
0.0000012300
0.0000015900
0.0000018700
-0.0000000241
0.0000090000
0.0000004540
0.0000000650
0.0000005440
0.0000001590
0.0000005530
0.0000002570
NaN
0.0000000354
0.0000012500
-0.0000000003
0.0000099700
0.0000000982
0.0000520000
0.0000000021
0.0000055000
0.0000082200
0.0000001180
0.0000009990
-0.0000002470
0.0000058000
0.0000196000
0.0000011200
NaN
0.0000023600
0.0000002270
0.0000022600
0.0000367000
0.0000011300
0.0000149000
0.0000134000
0.0000692000
0.0000282000
0.0000036500
0.0000033400
0.0000187000
0.0000000618
0.0000044700
0.0000028600
0.0001335000
0.0000005400
0.0000002470
0.0000104000
0.0000079500
0.0000054100
NaN
0.0000009290
0.0000006910
0.0000002660
-0.0000009680
0.0000011300
0.0000212000
0.0000009550
0.0000001540
0.0000006870
0.0000070600
0.0000027400
0.0000001090
0.0000001580
0.0000017800
0.0000006530
0.0000030400
0.0000017400
0.0000106000
0.0000048800
0.0000009890
0.0000043300
0.0000012000
0.0000001370
0.0000020100
0.0000000136
0.0000243000
0.0000015300
0.0000045500
0.0000035900
0.0000011100
0.0000003160
0.0000018400
0.0000087200
0.0000020700
0.0000037700
0.0000070300
0.0000040300
0.0000009220
0.0000017100
NaN
0.0000012600
0.0000016200
NaN
0.0000197000
0.0000025500
0.0000001280
0.0000366000
0.0000079600
0.0000013800
];

%% Industry markup and elasticity data
IMA=xlsread('Industry Markup Data.xlsx');
PME=NaN*ones(144,1);
for i=1:144
    IND=find(i==IMA(:,1))
    if isempty(IND)==0
        PME(i,1)=IMA(IND,19);
    end
end
TORQSTOCK(:,12)=PME;

%% Convert 1989-1990 CRSP stocks for use in comove tests
RAW=csvread(['ALLCRSPRETPRCVOL19891990REV.csv']); % load CRSP data on returns, prices, and volume over 1989-1990
PERMNOS=unique(RAW(:,2)); % find unique PERMNOs
P=length(PERMNOS); % number of PERMNOs
D1=find(FACTORS(:,1)==19891101); % Find first date index
D2=find(FACTORS(:,1)==19901031); % find last date index
TD=D2-D1+1; % find number of days
COMOVERETS=NaN*ones(TD,P); % Placeholder
COMOVEINDS=NaN*ones(P,1); % Placeholder
COMOVEMPRC=NaN*ones(P,1); % Placeholder
for p=1:P % loop over PERMNOs
    clc
    disp([num2str(p) ' of ' num2str(P)]) % on-screen counter
    IND=find(RAW(:,2)==PERMNOS(p)); % find the observations in CRSP data for each PERMNO
    if (length(IND)==253)&(sum(RAW(IND,7)<-1)==0); % if 253 days and no missing obs
        COMOVERETS(:,p)=RAW(IND,7); % Retrieve returns
        COMOVEINDS(p)=RAW(IND(end),3); % Retrieve industries
        COMOVEMPRC(p)=min(abs(RAW(IND,5)));  %minimum price over prior 253 days
    end
end
IND2=find(COMOVEINDS==0); % If industry is missing
COMOVERETS(:,IND2)=NaN; % replace with NaN
COMOVEINDS(IND2)=NaN;

%% Regress TORQ stocks on every other stock in CRSP over the year Nov 1, 1989 to Oct 31, 1990, using direct, Market Model, FF3F, FF4F controls.
TFIRMS=TORQSTOCK(:,2);
N=length(TFIRMS);
% OLS Coefficients
ATATRBR=NaN*zeros(N,P); % placeholder for r-bar
   
for t=1:N; % loop over sample firms (t)
    clear TORQID RETMATA
    clc
    disp(['Firm ' num2str(t) ' of ' num2str(N)]) % on screen counter
    TORQID=find(PERMNOS==TFIRMS(t)); % find sample firm in CRSP data
    RETMATA=COMOVERETS(:,TORQID); % pull return data for sample firm t
    for p=1:P; % loop over all PERMNOs in CRSP
        if p~=TORQID % do not include sample firm
            clear RETMATT REGDATA ATAT C CP
            RETMATT=COMOVERETS(:,p); % pull return data
            REGDATA=[RETMATA ones(253,1)  RETMATT];    %REGDATA is (250-Missing, 7) matrix that includes only good observations ZERODATE=(end-5), ESTDATE2=(end-11), EVDATE1=(end-10), EVDATE2=(end)
            ATAT=ols(REGDATA(:,1),REGDATA(:,[2 3])); %OLS Coefficients with Acq Rets = 1 + Tar Returns
            ATATRBR(t,p)=ATAT.rbar; % r-squared from regression
        end
    end
end

MAXRBRIND=NaN*ones(144,3); % placeholder
for t=1:length(TORQSTOCK(:,1)); % for each sample stock
    clear TEMP1 TEMP2
    TEMP1=flipdim(sortrows([ATATRBR(t,:)' COMOVEMPRC COMOVEINDS PERMNOS],1),1); % sort by r-squares
    IX1=find((TEMP1(:,1)<=1)&(TEMP1(:,2)>=5),1,'first'); % find maximum
    if isempty(IX1)==0
        MAXRBRIND(t,:)=TEMP1(IX1,[1 3 4]); % record maximum r2
    end
end
TORQSTOCK(:,13)=MAXRBRIND(:,1); % MAXRBR

%% Get risk-adjusted industry comovement
RAW=csvread(['ALLCRSPRET19891990REV.csv']);
PERMNOS=unique(RAW(:,3));
P=length(PERMNOS);
COMOVERETS=NaN*ones(TD,P);
COMOVEINDS=NaN*ones(P,1);
for p=1:P
   clc
   p
    IND=find(RAW(:,3)==PERMNOS(p));
    if (length(IND)==253)&(sum(RAW(IND,5)<-1)==0);
        COMOVERETS(:,p)=RAW(IND,5);
        COMOVEINDS(p)=RAW(IND(end),4);
    end
end
IND2=find(COMOVEINDS==0);
COMOVERETS(:,IND2)=NaN;
COMOVEINDS(IND2)=NaN;
load FACTORS
TFIRMS=TORQSTOCK(:,2);
N=length(TFIRMS);
P=length(PERMNOS);
%clear TORQSTOCK
COMOVEIND2=(COMOVEINDS-mod(COMOVEINDS,100))/100;
INDS2=unique(COMOVEIND2(COMOVEINDS<=9999));
F=length(COMOVERETS(1,:));
FF3FTORQCOGIND2MEDN=NaN*ones(N,1);
FF3FCOGIND2=NaN(length(COMOVERETS(1,:)),1);
FF3FCOGIND2MEDN=NaN(length(INDS2),1);
for f=1:length(COMOVERETS(1,:));
    clc
    f
    ID1=find((COMOVEIND2(f)==COMOVEIND2));
    if length(ID1)>5
        RETMATA=COMOVERETS(:,f);
        RETMATT=nanmean(COMOVERETS(:,ID1(ID1~=f)),2);
        REGDATA=[RETMATA ones(253,1) FACTORS(D1:D2,[2 3 6 7 8]) RETMATT]; %SAMMAT is the (250,8) [RETA 1 EW VW SMB HML UMD RETT] matrix of data for analysis
        FF3F=ols(REGDATA(:,1),REGDATA(:,[2 4 5 6 8])); %Fama French 3 Factor Model Coefficients + Tar
        FF3FCOGIND2(f,1)=FF3F.beta(5);
    end
    clear ID1 RETMATA RETMATT REGDATA FF3F
end

for i=1:length(INDS2);
    FF3FCOGIND2MEDN(i,1)=nanmedian(FF3FCOGIND2(find(INDS2(i)==COMOVEIND2)));
end

for t=1:length(TFIRMS);
    clear ID1 ;
    clc
    t
    ID1=find(PERMNOS==TFIRMS(t));
    if isnan(COMOVEINDS(ID1))==0
     FF3FTORQCOGIND2MEDN(t,1)=FF3FCOGIND2MEDN(find(INDS2==COMOVEIND2(ID1)));
    end
end
TORQSTOCK(:,14)=FF3FTORQCOGIND2MEDN;


%% Get Industry Comovement
ATATCOBIND2=NaN(F,1);
   
for f=1:length(COMOVERETS(1,:));
    clc
    disp([num2str(f) ' of ' num2str(length(COMOVERETS(1,:)))])
    ID1=find((COMOVEIND2(f)==COMOVEIND2));
    if length(ID1)>5
        RETMATA=COMOVERETS(:,f);
        RETMATT=nanmean(COMOVERETS(:,ID1(ID1~=f)),2);
        REGDATA=[RETMATA ones(253,1) RETMATT]; %SAMMAT is the (250,8) [RETA 1 EW VW SMB HML UMD RETT] matrix of data for analysis
        ATAT=ols(REGDATA(:,1),REGDATA(:,[2 3])); %Fama French 3 Factor Model Coefficients + Tar
        ATATCOBIND2(f,1)=ATAT.beta(2);
    end
    clear ID1 RETMATA RETMATT REGDATA ATAT

end

for i=1:length(INDS2);
    ATATCOBIND2MEDN(i,1)=nanmedian(ATATCOBIND2(find(INDS2(i)==COMOVEIND2)));
end

for t=1:length(TFIRMS);
    clear ID1 ;
    clc
    disp([num2str(t) ' of ' num2str(length(TFIRMS))])
    ID1=find(PERMNOS==TFIRMS(t));
    if isnan(COMOVEINDS(ID1))==0
         TORQSTOCK(t,15)=ATATCOBIND2MEDN(find(INDS2==COMOVEIND2(ID1)));
    end
end

%% This program collects data from IBES analyst forecasts to get COVA and add to TORQSTOCK.mat
IBES=xlsread(['IBES8990.xls']);
%TORQ Cusip Permno Match
TCP=xlsread(['TORQCUSIPPERMNO.xls']);
COVA=NaN*ones(144,1);
for i=1:144 % loop over sample firms
    clear CUSIP TDATA
    CUSIP=TCP(find(TORQSTOCK(i,2)==TCP(:,2)),1); % find CUSIP match
    TDATA=IBES(find(IBES(:,1)==CUSIP),:); %temporary data set
    COVA(i,1)=nanmean(TDATA(:,5)./abs(TDATA(:,4))); %Coefficient of variation of EPS (std dev/abs(mean))
end
TORQSTOCK(:,16)=COVA;
TORQSTOCK(find(TORQSTOCK(:,16)==Inf),16)=NaN; % replace infinity with missing

%% Compute R2 from FF model and saves per sample stock
for t=1:N;
    clear IDT RETMAT REGDATA LSVW FF3F FF4F
    IDT=find(PERMNOS==TFIRMS(t));
    RETMATA=COMOVERETS(:,IDT);
    REGDATA=[RETMATA-FACTORS(D1:D2,4) ones(253,1) FACTORS(D1:D2,3)-FACTORS(D1:D2,4) FACTORS(D1:D2,[6 7 8])]; %SAMMAT is the (250,8) [RETA-rf 1 VW-rf SMB HML UMD] matrix of data for analysis
    FF3F=ols(REGDATA(:,1),REGDATA(:,[2:5])); %Fama French 3 Factor Model Coefficients
    FF3FR2(t,1)=FF3F.rsqr;
end
TORQSTOCK(:,17)=FF3FR2;

%% Get Market Vigintiles
% These data are the vigintiles of each of the sample stocks based on all
% CRSP data.  I record them here, rather than including all CRSP/Compustat
% data for all stocks in the merged CRSP and Compustat database.
TORQSTOCK(:,18:21)=[2	7	NaN	NaN
1	5	14	19
NaN	NaN	NaN	NaN
1	1	1	1
17	19	19	15
NaN	NaN	NaN	NaN
5	11	6	3
13	10	10	6
15	12	1	15
NaN	NaN	NaN	NaN
18	20	14	12
20	1	1	10
16	16	13	14
18	17	4	5
6	5	20	18
11	8	11	10
6	18	2	10
18	20	NaN	NaN
NaN	NaN	NaN	NaN
NaN	19	4	1
2	9	19	15
1	2	9	13
12	20	9	10
11	18	9	19
12	11	9	16
3	5	15	10
19	1	2	7
18	12	19	10
2	7	19	18
3	7	11	13
4	1	10	19
13	1	5	4
1	19	10	2
18	20	12	10
8	12	13	12
12	10	11	16
14	14	14	8
7	6	17	16
17	1	1	1
1	3	10	15
NaN	NaN	NaN	NaN
NaN	NaN	1	1
20	1	NaN	NaN
19	16	3	13
19	1	7	20
13	20	13	6
9	9	16	9
NaN	18	14	4
20	1	2	14
NaN	NaN	NaN	NaN
5	15	20	8
7	15	10	8
17	15	17	15
17	14	18	13
17	1	2	2
5	9	20	16
NaN	NaN	1	1
1	6	NaN	NaN
12	11	16	6
5	8	5	7
14	1	2	6
1	13	NaN	NaN
NaN	NaN	NaN	NaN
13	11	12	11
14	15	7	20
17	19	3	6
NaN	NaN	NaN	NaN
12	9	20	18
20	1	NaN	NaN
13	15	9	17
15	19	NaN	NaN
19	1	1	4
20	1	1	1
NaN	NaN	NaN	NaN
1	1	14	11
12	18	5	2
9	20	NaN	NaN
12	18	14	6
14	12	9	14
NaN	NaN	NaN	NaN
9	4	NaN	NaN
9	18	4	8
13	14	12	16
9	5	11	6
NaN	NaN	NaN	NaN
3	11	20	17
NaN	NaN	3	1
18	18	NaN	NaN
4	5	NaN	NaN
8	13	10	7
18	15	15	14
8	1	4	4
NaN	NaN	NaN	NaN
7	9	15	14
NaN	NaN	1	1
15	13	17	12
8	9	19	17
11	13	6	6
13	1	1	4
8	1	13	18
NaN	NaN	NaN	NaN
11	12	15	14
NaN	NaN	NaN	NaN
3	7	7	2
1	5	10	17
13	11	11	4
14	14	17	14
16	17	18	17
NaN	NaN	10	20
NaN	NaN	17	20
16	1	4	10
18	19	NaN	NaN
3	1	NaN	NaN
14	12	8	12
2	3	20	20
7	1	1	1
11	6	17	20
NaN	NaN	NaN	NaN
NaN	NaN	14	20
3	6	20	8
20	1	NaN	NaN
16	3	11	12
15	17	NaN	NaN
19	13	7	8
9	6	8	10
15	15	18	16
15	1	2	3
1	1	1	1
18	16	11	7
7	10	12	13
11	6	8	7
6	1	1	12
3	1	9	20
15	17	16	11
NaN	NaN	1	2
NaN	NaN	8	2
5	6	11	7
20	1	NaN	NaN
1	20	1	1
3	5	17	19
10	4	20	18
16	14	12	12
NaN	NaN	NaN	NaN
13	17	8	9
];

%% Get Industry Maturity data
CS=csvread(['COMPUSTATSICS19501990.Coded.csv']); % Sales and SIC codes for all Compustat firms from 1950 to 1990, columns=GVKEY DATADATE FISCALYEAR SALES SIC
S=unique(CS(:,5));
SR=[sort(repmat(S,1990-1949,1)) repmat([1950:1990]',length(S),1) NaN*ones(length(S)*(1990-1949),1)];
for i=1:length(SR);
    clc
    disp([num2str(i) ' of ' num2str(length(SR))])
    ID=find((floor(CS(:,2)/10000)==SR(i,2))&(CS(:,5)==SR(i,1)));
    SR(i,3)=length(ID);
end

% find peak
for s=1:length(S);
    ID=find(SR(:,1)==S(s));
    [I J]=max(SR(ID,3));
    PEAK(s,:)=[S(s) SR(ID(J),2)];
end

MATURE=NaN*ones(144,1);
INDPEAK=NaN*ones(144,1);
for i=1:144;
    ID=find(SIC(i)==PEAK(:,1));
    if isempty(ID)==0
        MATURE(i,1)=(PEAK(ID,2)<1990)*1;
        INDPEAK(i,1)=PEAK(ID,2);
    end
end
TORQSTOCK(:,22:23)=[MATURE INDPEAK];

% Save a matlab file 
save(['TORQSTOCK'],'TORQSTOCK');

% Save a csv file to convert to Stata .dta file
STATA_HD={'symbol',...
'permno',...
'gvkey',...
'volatility',...
'turnover',...
'logassets',...
'elas_stock',...
'pin',...
'ffind5',...
'mrr',...
'lambda',...
'elas_prod',...
'maxr2',...
'riskadjcomove',...
'indcomove',...
'dispersion',...
'riskadjr2',...
'vigme',...
'vigpr',...
'vigbm',...
'vigep',...
'mature',...
'indpeak',...
};

cell2csv('TORQSTOCK.csv',STATA_HD,',');
dlmwrite('TORQSTOCK.csv',TORQSTOCK,'delimiter',',','precision','%1.9f','-append');

